%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Marginal Effect for model 4: with both effect, no unobs, without dummies%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% To study the marginal effect, consider an individual in a specific
%%%% group. So first, find a group

bnew1=bnew;
gg=71;  %% Find the groups whose mean characteristics is closest to the population.
            
X_f = dat(gg).xt;
n=size(X_f,1);
k=17;

bnew=bnew(1:2*k+2,:);
W_f= dat(gg).wt;

xx_new=dat2(gg).Mg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

indpv=[X_f W_f*X_f xx_new ones(n,1)];
F=1./(1+exp(-2*(indpv*bnew)));

naivemarg=[2*diag(F.*(1-F))*(bnew(1)+.2*indpv(:,1)*bnew(17)) 2*F.*(1-F)*[bnew(2:16)' bnew(18:35)']];
naivemargm=mean(naivemarg(1:n,:))'*100;  


    %%%%% Marginal effect for the discrete variables
ini_n=zeros(n,k-1);
    
X_fa=X_f; 
for member=1:n   

    for kk=3:k-1   %%%to k-1 ,as the last regressor is squared-age
        
        X_fa(member,kk)=0;
        indpv=[X_fa W_f*X_fa xx_new ones(n,1)];%
        F_before=1./(1+exp(-2*indpv*bnew));

        X_fa(member,kk)=1;

        indpv=[X_fa W_f*X_fa xx_new ones(n,1)];%
        F_after=1./(1+exp(-2*indpv*bnew));
       
        amarg=F_after-F_before;
        ini_n(member,kk)=amarg(member);

     
        X_fa=X_f;
    end
end
ini_n2=ini_n(:,3:k-1);
ma1_n=mean(ini_n2)'*100; %By multiply by 100, express in terms of percentage point
naivemargm2=naivemargm;
naivemargm2(3:k-1)=ma1_n;


fderiv_tanh=4./(exp(2*indpv*bnew)+exp(-2*indpv*bnew)+2); %%nx1

ini=zeros(n,k-1);
affect=zeros(n,k-1);
X_fa=X_f;  
for member=1:n 
    %%%%% Marginal effect for the continuous variables
    kk=2;  %%% years in school
    dXdx=zeros(size(X_f));
    dXdx(member,kk)=1;
    partial=inv(eye(n)-diag(fderiv_tanh)*bnew(2*k+1)*W_f)*diag(fderiv_tanh)*(dXdx*bnew(1:k)+W_f*dXdx*bnew(k+1:2*k)); %%% Calculate dM/db for group g
    amarg=2*diag(F.*(1-F))*(dXdx*bnew(1:k)+W_f*dXdx*bnew(k+1:2*k)+bnew(2*k+1)*W_f*partial); %% marginal effect for ALL members
    ini(member,kk)=amarg(member);
    amarg(member)=0;
    affect(:,kk)=affect(:,kk)+amarg;

    kk=1; %%% age
    dXdx=zeros(size(X_f));
    dXdx(member,kk)=1;
    dXdx(member,kk+16)=2/10*X_f(member,1);
    partial=inv(eye(n)-diag(fderiv_tanh)*bnew(2*k+1)*W_f)*diag(fderiv_tanh)*(dXdx*bnew(1:k)+W_f*dXdx*bnew(k+1:2*k)); %%% Calculate dM/db for group g
    amarg=2*diag(F.*(1-F))*(dXdx*bnew(1:k)+W_f*dXdx*bnew(k+1:2*k)+bnew(2*k+1)*W_f*partial); %% marginal effect for ALL members
    ini(member,kk)=amarg(member);
    amarg(member)=0;
    affect(:,kk)=affect(:,kk)+amarg;

   

    %%%%% Marginal effect for the discrete variables

    for kk=3:k-1   %%%to k-1 ,as the last regressor is squared-age

        X_fa(member,kk)=0;
        x0=ones(n,1)/2;
        euclid=1;
        while euclid>1e-7
            xx_new=tanh(X_fa*bnew(1:k)+W_f*X_fa*bnew(k+1:2*k)+W_f*x0*bnew(2*k+1)+bnew(2*k+2)*ones(n,1));
            euclid=max(sum((xx_new-x0).*(xx_new-x0)));
            x0=xx_new;
        end
        xx_new=W_f*xx_new;
        
        indpv=[X_fa W_f*X_fa xx_new ones(n,1)];%
        F_before=1./(1+exp(-2*indpv*bnew));

        X_fa(member,kk)=1;
        x0=ones(n,1)/2;
        euclid=1;
        while euclid>1e-7
            xx_new=tanh(X_fa*bnew(1:k)+W_f*X_fa*bnew(k+1:2*k)+W_f*x0*bnew(2*k+1)+bnew(2*k+2)*ones(n,1));
            euclid=max(sum((xx_new-x0).*(xx_new-x0)));
            x0=xx_new;
        end
        xx_new=W_f*xx_new;
        
        indpv=[X_fa W_f*X_fa xx_new ones(n,1)];%
        F_after=1./(1+exp(-2*indpv*bnew));
        clear xx_new
        clear euclid
        clear x0


        amarg=F_after-F_before;
        ini(member,kk)=amarg(member);
        amarg(member)=0;
        affect(:,kk)=affect(:,kk)+amarg;

   
        X_fa=X_f;
    end
end

ma1=mean(ini);  %%% average marginal effect of initiators
ma2=mean(affect)/n;%%% average marginal effect of followers
ma=[ma1; ma2]'*100 %%% By multiply by 100, express in terms of percentage point
[naivemargm naivemargm2];
naive=naivemargm2



